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Structure and bonding of dense liquid oxygen from first principles simulations 
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Using first principles simulations we have investigated the structural and bonding properties of 
dense fluid oxygen up to f 80 GPa. We have found that band gap closure occurs in the molecular 
liquid, with a "slow" transition from a semi-conducting to a poor metallic state occurring over a 
wide pressure range. At approximately 80 GPa, molecular dissociation is observed in the metallic 
fluid. Spin fluctuations play a key role in determining the electronic structure of the low pressure 
fluid, while they are suppressed at high pressure. 



The study of pressure-induced transformations in 
solids and liquids has become an active field of research 
in the last decade Q , due to key progress both in exper- 
imental methods (diamond anvil cell (DAC) and shock 
wave experiments) and in accurate computational tools. 
In particular, molecular systems such as hydrogen 0, 
oxygen |3( and nitrogen have been investigated with 
a variety of techniques. Nevertheless, fundamental ques- 
tions on the nature of these systems in e.g. the dense 
liquid state are yet unanswered. 

In this Letter, we focus on the liquid state of oxy- 
gen. In particular, we examine the structural and bond- 
ing changes induced by pressure, as well as the mecha- 
nism leading to the insulator-to-metal transition (ITMT) 
observed experimentally in shock wave experiments 0. 

Metallization of oxygen was first observed in the solid 
state by Desgreniers et al. who reported a Drude- 
type metallic behavior of its reflectivity for pressures 
(P) above 96 GPa. These results were later confirmed 
by Shimizu and coworkers who observed a change in 
the sign of the resistance-temperature curve at 96 GPa, 
and also discovered superconductivity in solid, dense 
oxygen ||. Following measurements by Akahama et 
al. 0, Week et al. 10] showed that in solid oxy- 
gen the ITMT is accompanied by a structural change 
between two different crystallographic arrangements of 
molecules. The insulating phase below 96 GPa was also 
studied with infrared absorption spectroscopy by Gorelli 
et al. 0] wno suggested the existence of a lattice of O4 
molecules at high pressure. However, this hypothesis is 
not confirmed by recent ab initio calculations |l2| , show- 
ing that a geometrical configuration built from extended 
"herringbone" -type chains is energetically favored over a 
lattice of O4 molecules. The structure predicted theoret- 
ically is consistent with earlier suggestions based on DAC 
experiments [l3| . 

In fluid oxygen, metallization was recently reported 
by Bastea et al. [5| who used reverberation shock wave 
techniques to generate compression states on a quasi- 
isentrope. For a given shock speed, the thermodynamic 
states reached in these experiments are at higher pres- 
sure and lower temperatures (T) than Hugoniot states 
obtained in earlier shock wave measurements [l4| . When 
increasing P from ~ 30 to 120 GPa (and correspondingly 



T by several thousands K), Bastea et al. [5| observed an 
increase of the electrical conductivity over six orders of 
magnitude. At about 120 GPa and 3.4-fold compression, 
the measured conductivity attains values characteristic 
of a poor metal. As P is further increased from 120 to 
190 GPa, the conductivity levels off and is nearly con- 
stant as a function of P. These experimental results, 
which indicate that the ITMT occurs at higher P in the 
liquid than in the solid, are rather surprising. In general, 
one would expect disorder or dissociation effects to lead 
to a lower metallization pressure in the fluid than in the 
ordered solid. 

In order to investigate the mechanism by which the 
ITMT occurs in the fluid, we have carried out a series 
of ab initio molecular dynamics (MD) simulations 0] in 
the liquid state for different pressures. Our findings shed 
light on the microscopic structure and electronic proper- 
ties of the dense fluid state and on the key role of spin 
fluctuations present in fluid before metallization occurs. 
Our predicted metallization pressure in the fluid is lower 
than that reported experimentally ||, as well as lower 
than the metallization pressure measured in the solid. 
Our findings further indicate that the system becomes 
metallic in a predominantly molecular state, as opposed 
to other fluids like hydrogen 0] and nitrogen [l7j . 

We have carried out ab initio MD simulations us- 
ing a gradient corrected energy functional 0] in the 
density (p) and temperature range explored in recent 
shock wave experiments: 2.8 < p < 4.5 gcm~ 3 and 
1000 < T < 6400 K. In some of our simulations we 
have treated the spin variable explicitly within the gradi- 
ent corrected local spin density approximation (GCSDA). 
We performed constant volume microcanonical simula- 
tions 0] with cells containing 54 and 108 oxygen atoms. 
For each p and T the total simulation time was 4-5 ps, 
after equilibration. 

Our calculated pair correlations functions {g(r)) are 
shown in Fig. ^ Finite size effects were found to be neg- 
ligible compared to temperature and pressure effects. At 
p = 2.8 gem" 3 and T = 1000 K [H, the first sharp peak 
of g(r) at 1.20 A, followed by a minimum very close to 
zero, is indicative of a molecular fluid. The computed 
pressure is 24 GPa. The first maximum is at a distance 
very close to the bond length of an oxygen molecule in the 
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FIG. 1: Pair correlation functions calculated using GCSDA 
(see text) for different densities, temperatures, and pressures 
as specified in the legends. 



gas phase in its triplet ground state (1.207 A) [2]]]. As the 
temperature is increased to 3000 K (and correspondingly 
the pressure to 34 GPa) the first molecular peak is broad- 
ened while its position remains unchanged. This trend 
continues as the fluid is compressed to p = 3.6 gem" 3 at 
T=3000K (P=84 GPa). Upon further compression to p 
= 4.5gcm" 3 at T=6400K (P=179 GPa), the first molec- 
ular peak becomes less intense and the second broad peak 
moves to shorter distances. It is interesting to note that 
the position of the intramolecular peak remains constant 
over a pressure range of 150 GPa and does not shift to 
larger separations as e.g. in dense liquid hydrogen [2*2^ . 

The ground state of the O2 molecule is a triplet, and 
this has been shown to give rise to interesting structural 
effects in the low density fluid (2^. It is therefore inter- 
esting to explore the effect of explicit spin treatment on 
the structural properties of the high pressure fluid as well. 
We found that spin effects are negligible on average struc- 
tural properties both at low and high pressure, consistent 
with recent ab initio calculations for solid oxygen |12J in 
the same pressure range. For example in Fig.^ g(r) cal- 
culated with and without explicit treatment of spin are 
found to be indistinguishable. However, spin effects are 
key in determining the optical excitation spectrum of the 
low-P fluid, as we will discuss below when we present our 
results for the electronic conductivity. 

The microscopic structure of the liquid can be further 
examined by a cluster analysis. We define a set of atoms 
as belonging to a given cluster if these atoms are sepa- 
rated by less than a given cut-off distance r c . (We chose 
r c = 1.50 A, which corresponds to the first g(r) minimum 
in the molecular regime.) We then define cluster lifetimes 
by simply computing the duration for which a given clus- 



ter remains intact. This analysis yields a very simple 
picture for our simulations at low density (2.8 gem -3 ) 
and 1000 and 3000 K: we find that the liquid consists en- 
tirely of O2 molecules. A density increase to 3.6 g cm" 3 at 
3000 K leads to a small fraction of dissociated molecules. 
Although the liquid is still predominantly molecular (94% 
O2), we find 1% single atoms and 5% 03-like clusters. At 
this density, the average lifetime of O2 molecules is about 
15 vibrational periods of the isolated O2 molecule (21.1 
fs), while Oi and O3 exist only for about one molecu- 
lar vibration. When the temperature is raised to 5000 K 
at 3.6 gem -3 , the number of molecules decreases to 79%, 
leading to 7% atoms, 10% O3, and 4% 04-like complexes. 
The lifetime of molecules decreases to about three O2 
vibrations, indicating the presence of a dynamical dis- 
sociation equilibrium with molecules still being the pre- 
dominant and most stable species. At p—4.5 gem" 3 and 
T=6400 K, we find 9% atoms and 58% molecules and also 
16% O3 clusters. The remaining 17% of the nuclei be- 
long to short-lived complexes composed of four or more 
atoms. Further inspection reveals that 04-like clusters 
arrange predominantly into short zig-zag chains, consis- 
tent with the finding of Ref. [12( for the compressed solid. 
At this high density, the O2 lifetime is only about 1.5 vi- 
brational periods. The short lifetime of O3 at all densities 
clearly indicates that triatomic configurations are closer 
to scattering states than to stable chemical species. Nev- 
ertheless, it is interesting to note that the "bond" angles 
of O3 clusters are not randomly distributed but vary be- 
tween 80° and 150° with an average of 120°. Therefore, 
the transient structure of O3 clusters resembles that of 
the ozone molecule. 




FIG. 2: Difference between spin-up and spin-down densities 
for two representative snapshots at p = 2.8 gem -3 and T = 
1000 K (left panel) and p = 4.5 gem" 3 and T = 6400 K (right 
panel). Positive and negative regions are represented by red 
and blue isosurfaces. Green spheres indicate the positions of 
the atoms. 



As mentioned above, while spin effects are negligible 
for the structural properties of the fluid, it is crucial to 
include them explicitly in order to understand the elec- 
tronic properties of liquid oxygen. In the low density 
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FIG. 3: Pair correlation function derived from nuclear posi- 
tions and the absolute difference of the spin density. The first 
peak quantitatively characterizes spin fluctuations present in 
the molecular liquid at low density of 2.8 gem" 3 . These fluc- 
tuations are suppressed at higher densities. 



regime, the liquid is composed of a disordered arrange- 
ment of O2 molecules carrying a net spin of one as in 
the gas phase. The left panel of Fig. |2l shows the differ- 
ence (Apjj) between spin-up (/Of) and spin-down densi- 
ties (pi), for a snapshot at p = 2.8gcm -3 and 1000 K, 
showing a random distribution of up and down spins. At 
higher P (right panel of Fig.EJ), A/?jj has almost vanished 
throughout the whole system, indicating that spin fluctu- 
ations are suppressed. The effect shown qualitatively in 
Fig. |21 can be characterized quantitatively by the correla- 
tion function of the nuclei and the electronic magnetic 
moment, m z (r) = p. B [pi{r) - Pi(r)}, 



Pego,mM r \) = J{ 



\m z {r l - r) 



(1) 



where p e is the electron density, N the number of atoms 
and (...) denotes an average over ionic configurations 
and the solid angle of r. This correlation function is sim- 
ilar to the nuclei-electron pair correlation function. The 
first, intense peak in the low density correlation function 
shown in Fig. [3] indicates that the majority of molecules 
are in the triplet state. At higher density, the intensity 
of the first peak is decreased by almost an order of mag- 
nitude, showing that most of the pairs in the fluid do not 
carry a net spin. 

In order to qualitatively understand the effect of pres- 
sure on the spin state of the oxygen molecule, we carried 
out calculations of the total energies of O2 in the sin- 
glet (Es=o) and triplet (Es=i) state in the presence of a 
confining potential mimicking an average pressure effect. 
Interestingly, we found that the difference Es=o — E$=i 
can be substantially reduced by applying an external con- 
fining potential. This difference decreases as a function 
of the confinement strength. 

In order to compare our results with recent shock wave 
data, we computed the conductivity (a) of the liquid for 



selected snapshots taken from our GCSDA simulations 
at differe nt p ressures, using the Kubo-Greenwood for- 
mula |3; l25| . We note that in the low pressure regime, 
where spin fluctuations are present, GCSDA and GGA 
results are very different. In particular, the GGA yields 
much smaller optical gaps, consistent with the incorrect 
treatment of excitation in the isolated molecule if spin is 
not explicitly taken into account. On the other hand, at 
higher densities, the computed GGA and GCSDA con- 
ductivities are very similar. 
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FIG. 4: Measured and calculated conductivity as a function of 
pressure. The labels on the simulation data points correspond 
to the thermodynamic conditions specified in Fig. 

Our results are reported in Fig. 0] along with the ex- 
perimental results of Ref. H and 0] . The comparison 
between theory and experiment is not straightforward: 
the temperature in the measurements by of Hamilton et 
al. [ljl is in general different from that of Bastea et al.'s 
data |5j, and in both experiments, T was not measured 
but estimated on the basis of empirical equation of state 
(EOS) models. However, the two sets of experimental 
data show the same qualitative rise of a as a function of P 
and T. In particular, Bastea et aZ.'s data show a six order 
of magnitude rise in conductivity with increasing shock 
compression, corresponding to a P and T increase from 
about 30 to 120 GPa and 1000 to 5000 K, respectively. 
Our calculations within GCSDA show only two orders 
of magnitude increase in a in the same P and T range. 
This increase of a as a function of T occurs primarily in 
a state where the fluid is entirely composed of stable O2 
molecules (p = 2.8gcm- 3 ,T = 1000 -> 3000 K). This 
indicates that in oxygen, metallization is not accompa- 
nied by molecular dissociation as found in other first row 
molecular fluids, in particular hydrogen |25j and nitro- 
gen ^t|- Above 100 GPa both theory and experiment 
show an almost constant a as a function of compression, 
but the theoretical values are about a factor of 7 larger 
than the experimental ones. 
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In our calculations, overestimates of conductivity val- 
ues could be due to the use of the GCSDA, which is 
known to underestimate electronic gaps in many systems. 
In order to qualitatively investigate this effect, we used 
a scissors operator and arbitrarily increased the value of 
the computed gaps by a factor of two. This lowered the 
computed conductivity by one order of magnitude at p = 
2.8 gem" 3 , both at 1000 and 3000 K. However, even after 
this "gap correction" , our theoretical data show values of 
a approaching those of a poor metal at a pressure much 
lower (30-40 GPa) than found experimentally. The com- 
puted conductivity is then slowly changing as a function 
of P, and it reaches metallic-like values at approximately 
80 GPa. This "slow" metallization of the fluid as a func- 
tion of pressure may be related to spin fluctuations and 
their "slow" suppression as a function of pressure. At 
low P, the molecular fluid, with each oxygen molecule in 
a triplet state, can be viewed as a disordered spin system. 
As the pressure is increased, this system shows a glass- 
like behavior due to the presence of a multitude of local 
minima in its potential energy surface, corresponding to a 
multitude of different spin arrangements. Eventually, the 
fluid reaches metallization when spin fluctuations are al- 
most completely suppressed. In this respect, compressed 
liquid oxygen is unique among first row molecular fluids 
and rather different from, e.g. liquid hydrogen and liquid 
nitrogen, composed of molecules with spin zero. 

In summary, our first principles simulations show that 
above 30 GPa, over a wide pressure range of about 50 
GPa, liquid oxygen transforms from a semi-conducting to 
a metallic fluid by closure of the band gap in a disordered 
molecular state. For P > 80 GPa, molecular dissociation 
is observed, giving rise to short-lived atomic- and ozone- 
like species. We did not observe any stable O4 molecule 
for the thermodynamic conditions considered here, which 
is consistent with recent theoretical findings for solid oxy- 
gen 01 • O ur results also show that spin fluctuations are 
suppressed at high pressure, with no significant magnetic 
moment present in the high pressure liquid, as opposed 
to the semi-conducting low pressure regime. 

Our findings for the electrical conductivity point at 
the need for further experimental and theoretical work. 
Experimentally, in shock wave experiments only P and 
a are measured. Data are then analyzed using a model 
EOS as input for hydrodynamic codes which compute 
densities. Finally, the temperature is inferred from addi- 
tional empirical chemical model calculations. This type 
of analysis may introduce uncertainties in the estimated 
p and T at a given P and a. Measurements of the oxy- 
gen EOS obtained with shock wave experiments and pos- 
sibly with other techniques, e.g. DAC measurements, 
may help clarify the existing discrepancies between the- 
ory and experiment. Finally, we note that theoretical 
inaccuracies on calculated conductivity values may arise 
from the use of an approximate exchange and correlation 
functional. Although computationally still out of reach, 



quantum Monte Carlo calculations for liquid oxygen un- 
der pressure might be used in the future to address this 
issue (lif . 
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